Plantilla supervivència
VHIR-UEB-for-019v.01
1 Introducció
Aquesta plantilla mostra un exemple d’anàlisi de supervivència, tècnica estadística utilitzada per a analitzar el temps que triga a ocórrer un esdeveniment, com la mort, la falla d’un component o l’èxit en un tractament mèdic. És a dir, l’anàlisi de supervivència avalua i prediu el temps fins que es produeix un succés.
Els objectius són: calcular la probabilitat que no succeeixi l’esdeveniment en un temps T, comparar la supervivència entre grups i analitzar factors de risc associats a la supervivència.
Per a poder realitzar aquest anàlisi hem de tenir la següent informació per a cadascun dels individus: l’esdeveniment que es vol analitzar; una data origen, que pot ser la data d’inclusió o la data de diagnòstic; una data final, que pot ser la data que ha ocorregut l’esdeveniment o la última data disponible; el període de temps fins que ha aparegut aquest esdeveniment; i informació addicional.
A vegades, al final del seguiment, alguns dels individus no han tingut l’esdeveniment d’interès i, per tant, es desconeix el seu veritable temps fins a l’esdeveniment, això és el que es coneix com a censura.
Ens hem d’assegurar que les variables estiguin definides de la següent manera: esdeveniment, amb valor 0 si l’individu ha estat censurat, s’ha recuperat o no es coneix el resultat final, i valor 1 si s’ha observat l’esdeveniment; i temps de seguiment, com el temps entre el moment d’inici i l’aparició de l’esdeveniment o final del seguiment.
Abans de començar amb aquesta anàlisi en programari R, es recomana consultar la plantilla de descriptius on estan explicats els passos previs de qualsevol anàlisi. A continuació, ja es pot procedir a l’anàlisi de supervivència.
En aquest exemple utilitzarem la base de dades diabetes.
## Warning: package 'knitr' was built under R version 4.4.1Tots els anàlisis han estat realitzats amb el programa estadístic “R”( R version 4.4.0 (2024-04-24 ucrt), Copyright (C) 2015 The R Foundation for Statistical Computing ).
> ##################
> #### GRUPS DE VARIABLES
> ##################
>
> var_id <- "id"
> var_event <- "mort"
> data_event <- "date_event"
> data_ini <- "date_ini"
> var_time <- "tempsviu"
> varsurv <- names(dat %>% select(edat,bmi,edatdiag,tabac,sbp,dbp,ecg,chd))
>
> ##################2 Temps a mort
En aquest apartat es mostra, de forma interactiva, les dades de supervivència individuals tenint en compte les dates reals de seguiment. Aquest primer apartat és molt importat per a poder visualitzar els temps fins a la mort de manera general.
La funció ggplotly() del paquet plotly a partir d’un ggplot() crea un gràfic interactiu el qual una vegada creat pots descarregar-ho en format png, apropar o allunyar a una zona concreta, moure’t pel gràfic, etc, a més mostra per a cada individu una línia corresponent al temps, l’esdeveniment o punt final i una caixa de text amb: Identificador, data inici de seguiment, data de l’esdeveniment o últim seguiment i estatus vital.
> ggplotly(
+ dat %>%
+ mutate(text = paste("ID pacient= ", !!sym(var_id), "<br>",
+ "Data inici = ", !!sym(data_ini), "<br>",
+ "Data event = ", !!sym(data_event), "<br>",
+ "Edat = ", edat)) %>%
+ ggplot(aes_string(x = "id", y = data_event, text = "text")) +
+ geom_linerange(aes_string(ymin = data_ini, ymax = data_event)) +
+ geom_point(aes_string(shape = var_event, color = var_event), stroke = 1, cex = 2) +
+ scale_shape_manual(values = c(1, 3, 4)) +
+ scale_color_manual(values = c("#20B2AA","#8A2BE2")) +
+ labs(y = "Temps", x = "ID pacient") +
+ coord_flip() +
+ theme_bw(base_size = 8) ,
+ tooltip = "text"
+ )Aquest mateix gràfic es realitza ajustant els valors d’inici de seguiment a 0, així es poden comparar millor els temps. Per a cada individu s’indiquen diverses característiques: Identificador, temps de seguiment (en anys), estat i edat.
> ggplotly(
+ dat %>%
+ mutate(
+ text = paste("ID pacient= ", !!sym(var_id),"<br>",
+ "Temps = ", round(!!sym(var_time),2), "<br>",
+ "Event = ", !!sym(var_event), "<br>",
+ "Edat = ", round(edat,1)) ) %>%
+ ggplot(aes_string(x = "id", y = "tempsviu", text = "text")) +
+ geom_linerange(aes_string(ymin = 0, ymax = "tempsviu")) +
+ geom_point(aes_string(shape = "mort", color = "mort"), stroke = 1, cex = 2) +
+ scale_shape_manual(values = c(1, 3, 4)) +
+ scale_color_manual(values = c("#20B2AA","#8A2BE2")) +
+ labs(y = "Temps (Anys)", x = "ID pacient") +
+ coord_flip() +
+ theme_bw(base_size = 8),
+ tooltip = "text")3 Anàlisi de supervivència
En aquest apartat s’ha realitzat un anàlisi de supervivència global i l’estimació de la corba Kaplan-Meier.
Primer es crea un objecte Surv amb la funció Surv(Time,Event) del paquet survival, que pren com a arguments: Time com al temps fins a l’esdeveniment i Event com a indicador si s’ha produït l’esdeveniment. El resultat és un vector del temps de seguiment, amb el símbol “+” per a representar aquelles observacions censurades.
La funció de supervivència ens indica la probabilitat de sobreviure (no ocorre l’esdeveniment) més enllà d’un període de temps t: \(S(t) = Prob(\text{Sobreviure t}) = P\{T>t\} = \int^{\infty}_{t}f(t) du= 1-F(t)\)
On: - Funció de densitat: \(f(t) = Prob(\text{instantània de morir})\) - Funció de distribució: \(F(t) = Prob(\text{morir abans de t})\)
A continuació, la funció surfit() del paquet survival retorna l’estimació de la corba de supervivència global, per veure un resum utilitzem summary() que retorna els següents valors: time, temps d’observació; n.risk, el nombre d’individus en risc; n.event, el nombre d’individus que presenten l’esdeveniment; survival, l’estimació de la funció de supervivència; std.err, la desviació estàndard de l’estimació; i lower/upper 95% CI, els intervals de confiança per a l’estimació. Per a l’estimació global la fórmula ha de ser Surv(time,event) ~ 1, per defecte type = “kaplan-meier, però l’estimació també pot ser d’altres tipus.
L’estimació no paramètrica (distribució del temps de supervivència desconeguda) de Kaplan-Meier de S(t) que assumeix que l’esdeveniment és independent per a cada individu, es calcula de la següent manera: \(\hat{S}(t)=\prod_{t_i<t}\frac{n_i - d_i}{n_i}\)
On:
- \(d_i\), número de morts en el moment t_i
- \(n_i\), número d’individus en risc abans de t_i. Si no hi ha censura és el número de supervivents immediats abans del mont t_i, amb censura és el número de supervivents menys el número de casos censurats.
Per últim és grafica la corba de supervivència global estimada amb la funció ggsurvplot() del paquet survminer. La corba de supervivència pot ser transformada amb fun = c(“event”, “cumhaz”, “pct”), per a esdeveniments acumulats, risc acumulat i corba de supervivència en percentatge.
> # guardem l'objecte Surv
> my.surv <- Surv(time = dat[,var_time], event = as.numeric(dat[,var_event]) - 1)
> head(my.surv, 10)## [1] 12.4+ 12.4+ 9.6+ 7.2+ 14.1+ 14.1+ 12.4+ 14.2+ 12.4+ 14.5+> # estimació Kaplan Meier
> fit_surv_glob <- survfit(my.surv ~ 1, data = dat )
>
> (tt <- t(data.frame(SurvGlobal = summary(fit_surv_glob)$table[!names(summary(fit_surv_glob)$table) %in% c("records","n.max","*rmean","*se(rmean)")])))## n.start events rmean se(rmean) median 0.95LCL 0.95UCL
## SurvGlobal 149 25 14.99397 0.3464687 NA NA NA> #xtable(tt, caption = "label{smuerte}Análisis de supervivència " )
>
> # estimació de la funció de supervivència
> summary(fit_surv_glob)## Call: survfit(formula = my.surv ~ 1, data = dat)
##
## time n.risk n.event survival std.err lower 95% CI upper 95% CI
## 2.0 142 1 0.993 0.00702 0.979 1.000
## 2.2 141 1 0.986 0.00989 0.967 1.000
## 3.6 137 2 0.972 0.01404 0.944 0.999
## 4.8 130 1 0.964 0.01579 0.934 0.996
## 5.5 128 2 0.949 0.01880 0.913 0.987
## 6.5 123 1 0.941 0.02017 0.903 0.982
## 6.7 122 1 0.934 0.02143 0.892 0.977
## 7.0 118 1 0.926 0.02266 0.882 0.971
## 7.2 117 1 0.918 0.02381 0.872 0.966
## 7.6 111 1 0.909 0.02499 0.862 0.960
## 8.0 110 1 0.901 0.02609 0.851 0.954
## 9.8 103 1 0.892 0.02727 0.841 0.948
## 10.0 102 1 0.884 0.02837 0.830 0.941
## 10.6 92 1 0.874 0.02964 0.818 0.934
## 10.8 89 1 0.864 0.03089 0.806 0.927
## 10.9 87 1 0.854 0.03209 0.794 0.920
## 11.0 85 1 0.844 0.03325 0.782 0.912
## 12.1 70 1 0.832 0.03490 0.767 0.904
## 12.5 52 1 0.816 0.03772 0.746 0.894
## 13.6 43 1 0.797 0.04134 0.720 0.883
## 13.7 41 1 0.778 0.04467 0.695 0.870
## 14.3 28 1 0.750 0.05099 0.656 0.857
## 15.3 13 1 0.692 0.07272 0.564 0.851> # grafic de la corba de KM
> ggsurvplot(fit_surv_glob, data = dat, risk.table = T, xlab = "anys")
> #plot(fit_surv_glob, mark.time = T, xlab = "anys", main = paste("Kaplan-Meier estimate with 95% confidence bounds.", var_event), cex.main = 0.8 )
>
> #(tt <- t(data.frame(SurvGlobal = summary(fit_surv_glob)$table[!names(summary(fit_surv_glob)$table) %in% c("records","n.max","*rmean","*se(rmean)")])))
> # xtable(tt, caption = "label{smuerte}Análisis de supervivència " )Una forma fàcil d’interpretar el gràfic és dir que en el moment zero, tots els participants segueixen vius i la probabilitat de supervivència és del 100%. Aquesta probabilitat disminueix amb el temps a mesura que moren els pacients. Els símbols “+” representen les observacions censurades o que han sortit de l’estudi.
La mediana de supervivència es pot estimar amb el percentil 50 de la distribució, que correspon amb el primer temps amb una proporció de supervivència igual o inferior a 0.5. Amb el gràfic només fa falta traçar una línia horitzontal a 0.5 i veure a quants anys fa referència. Com en aquest exemple la probabilitat de sobreviure no baixa del 60%, podem dir que la proporció de participants que no moren més enllà dels 15 anys de seguiment se situa més o menys al 70%.
3.1 Anàlisi univariant
> pval_cut <- 0.05Per a les variables qualitatives s’ha estimat la corba de supervivència Kaplan-Meier per a cadascuna de les categories per a saber si hi ha diferències entre elles. Per a l’estimació la fórmula ha de ser Surv(time,event) ~ grup.
En cada gràfic es mostren les diferents corbes juntament amb el valor p resultant de la prova log-rank o test de riscos proporcionals de comparació de corbes de supervivència. Si el p \(<\) 0.05 indica que hi ha diferències en la supervivència (si es compleixen les condicions de proporcionalitat de riscos).
Aquesta prova estadística de contrast d’hipòtesis compara el número d’esdeveniments en cada grup amb el número d’esdeveniments esperats combinant els diferents grups i és útil per detectar diferències a llarg termini.
La fórmula de l’estadístic de Log-rank per comparar dos grups és: \(\text{Log-Rank (aproximació)} = \frac{(O-E)^2}{(Var(O-E))}\), on O són els esdeveniments observats i E els esperats.
Si es comparen més de dos grups la fórmula es: \(\text{Log-Rank (aproximació)} = \sum^k_i\frac{(O_i-E_i)^2}{E_i}\), on k és el número de grups.
Contrast: \[ \left\{ \begin{array}{ll} H_{0}: & \text{les supervivències dels grups son la mateixa} \\ H_{1}: & \text{al menys un dels grups té una supervivència diferent} \end{array} \right. \]
A més es mostra una taula resum que inclou el nombre de casos, el nombre d’esdeveniments observats, la mediana en anys i els seus intervals de confiança al 95%.
Per a les variables que presenten diferències entre grups, s’ha ajustat un model de Cox univariant amb la funció coxph() del paquet survival, que s’utilitza quan es vol estimar l’efecte d’un conjunt de factors en la supervivència. La taula resultant indica el Hazard Ratio i el corresponent interval de confiança junt amb al valor p resultant de la prova log-rank.
El Hazard Ratio és el quocient entre els riscos dels diferents grups, i expressa la probabilitat que un individu experimenti un esdeveniment en un determinat moment, assumint que aquest individu ha sobreviscut fins a aquest instant sense experimentar l’esdeveniment d’objecte d’estudi. Un hazard ratio (HR) > 1 significa que l’exposició al factor augmenta la velocitat de produir-se l’esdeveniment, HR < 1 disminueix la velocitat, i HR = 1 diem que el factor no influeix en la supervivència.
La fórmula del Hazard Ratio es: \(\lambda(t)=h(t) = Prob(\frac{\text{morir entre t,t}+\Delta t}{\text{sobreviure a t}})=f(t)/S(t)\), i s’expressa amb un interval de confiança que és el rang de valors que probablement inclou el valor poblacional real.
> var2surv_cat <- c(names(which(unlist(lapply(dat[,c(varsurv)], is.factor)))))
> #var2surv_cat <- c(var2surv_cat[!var2surv_cat %in% "mort"])
> var2surv_num <- names(which(!unlist(lapply(dat[,c(varsurv)], is.factor))))
>
>
> splots <- list()
> fits <- list()
> pvalues <- NULL
>
> for (i in seq_along(var2surv_cat)) {
+
+ nm <- ifelse(Hmisc::label(dat[,var2surv_cat[i]]) != "",
+ Hmisc::label(dat[,var2surv_cat[i]]), var2surv_cat[i] )
+ cat(" \n###",nm, " \n")
+
+ #estimació kaplan-meier
+ fits[[i]] <- survfit(as.formula(paste("my.surv ~ ", var2surv_cat[i])), data = dat )
+ namevar <- ifelse(label(dat)[names(dat) == var2surv_cat[i]] != "",
+ label(dat)[names(dat) == var2surv_cat[i]], var2surv_cat[i])
+ # grafic amb p-valor
+ splots[[i]] <- ggsurvplot(fits[[i]], data = dat, pval = TRUE, risk.table = T,
+ xlab = "anys", title = paste0("Kaplan-Meier estimate (",namevar, ")."),
+ font.main = 13, risk.table.col = "strata", conf.int = TRUE,
+ legend.labs = levels(dat[,var2surv_cat[i]]))
+ print(splots[[i]])
+ print(xtable(summary(fits[[i]])$table[,c("records", "events", "median", "0.95LCL",
+ "0.95UCL")]), size = "small", type = "html")
+ pvalues[var2surv_cat[i]] <- surv_pvalue(fits[[i]])$pval
+
+ if (surv_pvalue(fits[[i]])$pval < pval_cut) {
+ ## fit a Cox model
+ mod <- coxph(as.formula(paste('my.surv~',var2surv_cat[i])), data = dat)
+ sum.surv <- summary(mod)
+ c_index <- sum.surv$concordance
+ ## Make pretty summary
+ print(kable(prettify(summary(mod), digits = 4)[,c(" ", "Hazard Ratio", "CI (lower)",
+ "CI (upper)","Pr(>|z|)" )], booktabs = T,
+ caption = paste0( "Hazard Ratio Cox ", namevar ,". Surv"),
+ longtable = TRUE, escape = F,col.names = c(" ", "Hazard Ratio", "CI (lower)",
+ "CI (upper)","P-value" )) %>%
+ kable_styling(latex_options = c("striped","hold_position", "repeat_header"),
+ font_size = 12) %>%
+ row_spec(0,background = "#993489", color = "white") )
+ print(ggforest(mod))
+ }
+ cat(" \n")
+ # print(var2surv_cat[i])
+ }3.1.1 tabac
| records | events | median | 0.95LCL | 0.95UCL | |
|---|---|---|---|---|---|
| tabac=No fumador | 57.00 | 11.00 | 15.30 | ||
| tabac=Ex fumador | 41.00 | 9.00 | |||
| tabac=Fumador | 51.00 | 5.00 |
3.1.2 ecg
| records | events | median | 0.95LCL | 0.95UCL | |
|---|---|---|---|---|---|
| ecg=Normal | 111.00 | 14.00 | |||
| ecg=Anormal | 11.00 | 7.00 | 7.20 | 5.50 | |
| ecg=Frontera | 27.00 | 4.00 |
| Hazard Ratio | CI (lower) | CI (upper) | P-value | |
|---|---|---|---|---|
| ecg: Anormal | 31.52 | 10.15 | 97.81 | <0.001 |
| ecg: Frontera | 2.356 | 0.7462 | 7.439 | 0.144 |
3.1.3 chd
| records | events | median | 0.95LCL | 0.95UCL | |
|---|---|---|---|---|---|
| chd=No | 99.00 | 13.00 | |||
| chd=Si | 50.00 | 12.00 | 12.50 |
| Hazard Ratio | CI (lower) | CI (upper) | P-value | |
|---|---|---|---|---|
| chd: Si | 3.542 | 1.559 | 8.045 | 0.003 |
En aquest exemple veiem que entre els grups de la variable tabac no hi ha diferències significatives en les corbes amb un p-valor de 0.15, en canvi, per a la variable ecg (lectura d’electrocardiograma) i chd (malaltia coronaria), si hi ha diferències, amb p-valors menors a 0.05. Per tant, podem dir que els que tenen una ecg normal tenen una probabilitat de sobreviure major que els que tenen una ecg anormal, i que els que no tenen chd tenen més probabilitat de sobreviure que els que sí que en tenen chd.
Per a les variables quantitatives s’ajusta un model de Cox univariant. La taula i el gràfic indiquen el Hazard Ratio i el seu corresponent interval de confiança. En el peu del grafic es pot observar el nombre total d’esdeveniments, el p-valor resultant de la prova log-rank i l’índex de concordança (probabilitat que, per a un parell de pacients comparables escollits a l’atzar, el pacient amb probabilitat de risc superior experimenti un esdeveniment abans que el pacient amb risc inferior).
> for (i in seq_along(var2surv_num)) {
+ namevar <- ifelse(label(dat)[names(dat) == var2surv_num[i]] != "",
+ label(dat)[names(dat) == var2surv_num[i]], var2surv_num[i])
+ ## fit a Cox model
+ mod <- coxph(as.formula(paste('my.surv~',var2surv_num[i])), data = dat)
+ sum.surv <- summary(mod)
+ c_index <- sum.surv$concordance
+ ## Make pretty summary
+ pvalues[var2surv_num[i]] <- summary(mod)$waldtest["pvalue"]
+
+ print(kable(prettify(summary(mod), digits = 4)[,c(" ", "Hazard Ratio", "CI (lower)",
+ "CI (upper)","Pr(>|z|)" )], booktabs = T,
+ caption = paste0( "Hazard Ratio Cox ", namevar ,". Surv"),
+ longtable = TRUE, escape = F, col.names = c(" ", "Hazard Ratio", "CI (lower)",
+ "CI (upper)","P-value" )) %>%
+ kable_styling(latex_options = c("striped","hold_position", "repeat_header"),
+ font_size = 12) %>%
+ row_spec(0,background = "#993489", color = "white") )
+ print(ggforest(mod))
+ }| Hazard Ratio | CI (lower) | CI (upper) | P-value | |
|---|---|---|---|---|
| edat | 1.102 | 1.067 | 1.139 | < 0.001 |
| Hazard Ratio | CI (lower) | CI (upper) | P-value | |
|---|---|---|---|---|
| bmi | 0.9482 | 0.8866 | 1.014 | 0.121 |
| Hazard Ratio | CI (lower) | CI (upper) | P-value | |
|---|---|---|---|---|
| edatdiag | 1.080 | 1.047 | 1.114 | < 0.001 |
| Hazard Ratio | CI (lower) | CI (upper) | P-value | |
|---|---|---|---|---|
| sbp | 1.024 | 1.005 | 1.045 | 0.015 |
| Hazard Ratio | CI (lower) | CI (upper) | P-value | |
|---|---|---|---|---|
| dbp | 1.022 | 0.9854 | 1.060 | 0.241 |
En aquest exemple les variables edat i edat del diagnòstic tenen p-valors significatius, per tant, a l’augmentar l’edat, el risc de morir creix un 10% aproximadament. En canvi, per a les variables bmi, sbp i dbp, no podem dir que les diferències siguin significatives.
3.2 Anàlisi multivariant
En aquest apartat es valorarà de manera simultània l’efecte d’una sèrie de variables explicatives o factors pronòstics en la supervivència.
3.2.1 Criteri estadístic
> complete_obs <- 90
> library(glmnet)
>
> vardel <- NAperc(dat[,varsurv],maxNA = 100 - complete_obs)$var
>
> var_mod_n <- varsurv[!varsurv %in% vardel]El següent model s’ha ajustat per criteri purament estadístic. S’han tingut en compte totes las variables que es mostren en l’ajust univariant i que el 90 % dels pacients completa. Les variables que s’han avaluat són: edat, bmi, edatdiag, tabac, sbp, dbp, ecg, chd
Per a seleccionar el model s’ha utilitzat la técnica Lasso de selecció de variables. És una tècnica que a l’hora d’ajustar penalitza la funció de maximització dels coeficients de les variables a partir d’un model on totes les variables s’ajusten fins a un final on tots els coeficients són 0. Aquesta penalització es regeix per un paràmetre lambda que es maximitza fins que es troba el valor òptim i les variables a ajustar. Aquesta tècnica permet controlar el sobreajust (overfitting) de forma que no s’inclouen variables supèrflues.
Par a efectuar-la se segueixen tres passos. En primer lloc, apliquem un ajust lasso par a veure que variables entrar i el número òptim. Després validem aquest procés amb validació creuada, és a dir, dividim la mostra en 10 parts de les quals escollim 9 com a entrenament i una com a validació i efectuem el procés anterior de Lasso 10 vegades i veiem les variables que millor expliquen la variable resultat. En tercer lloc, utilitzem un algoritme corregit del Lasso més eficient per a comprovar de nou que les variables finalment seleccionades són òptimes.
> # eliminem els NA
> d <- na.omit(dat[,c(var_mod_n, var_event, var_time)])
> # creem variable 01 da la variable resposta
> d[,paste0(var_event,"01")] <- as.numeric(d[,var_event]) - 1
> # canviem els temps 0 per 0.00001
> d[,var_time][which(d[,var_time] == 0)] <- 0.00001
>
> # eliminem el temps i les variables esdeveniment de les variables explicatives
> frml <- as.formula(paste0(" ~ . -", var_time, "- ", var_event, "-",paste0(var_event,"01")))
> #definim la matriu del model on les variables qualitatives estan transformades a numeriques
> x <- model.matrix( frml, d)
> # definim la supervivència
> y <- Surv(time = d[,var_time],
+ event = d[,paste0(var_event,"01")])
>
> # model amb regularització lasso
> fit <- glmnet(x, y, family = "cox", alpha = 1)
> plot(fit, label=T)
> # validació creuada
> set.seed(1)
> cv.fit <- cv.glmnet(x, y, family = "cox", alpha = 1)
> #seleccionem la lamda amb error minim
> plot_glmnet(fit,label = TRUE, s = cv.fit$lambda.min)
> #agamen els coeficients per a la lambda optima
> lasso_coef <- coef(fit, s = cv.fit$lambda.min)
> lasso_coef11 x 1 sparse Matrix of class “dgCMatrix”
1
(Intercept) .
edat 0.063219151
bmi .
edatdiag .
tabacEx fumador .
tabacFumador .
sbp .
dbp 0.004837547
ecgAnormal 1.879901596
ecgFrontera .
chdSi .
> #extraiem les variables seleccionades
> lasso_coef0 <- rownames(lasso_coef)[lasso_coef[,1] != 0]
> (names_lasso0 <- unique(extract_names(varlev = lasso_coef0,var = varsurv)))[1] “edat” “dbp” “ecg”
La regressió de Cox és una tècnica multivariant que permet identificar i avaluar la relació entre un conjunt de variables explicatives. També permet predir les probabilitats de supervivència per a un determinat individu a partir dels valors que prenen les seves variables pronostiques.
L’equació del model de regressió de Cox és: \(Ln(\lambda_t)=\beta_0+\beta_1x_1+\beta_2x_2+...+\beta_rx_r\)
Una vegada tenim les variables seleccionades realitzem el mateix model que en l’anàlisi univariant.
> frml <- formula(paste0("my.surv ~ ", paste0(names_lasso0, collapse = " + ") ))
> # model de Cox amb les variables seleccionades
> mod <- coxph(frml ,data = dat)
> sum.surv <- summary(mod)
> c_index <- sum.surv$concordance
>
> pvalues <- summary(mod)$waldtest["pvalue"]
>
> print(kable(prettify(summary(mod), digits = 4)[,c(" ", "Hazard Ratio", "CI (lower)", "CI (upper)","Pr(>|z|)" )], booktabs = T,
+ caption = paste0( "Hazard Ratio Cox Multivariate" ,". Surv. C Index:",round(c_index["C"],2)), longtable = TRUE, escape = F, col.names = c(" ", "Hazard Ratio", "CI (lower)", "CI (upper)","P-value" )) %>%
+ kable_styling(latex_options = c("striped","hold_position", "repeat_header"), font_size = 12) %>%
+ row_spec(0,background = "#993489", color = "white") )| Hazard Ratio | CI (lower) | CI (upper) | P-value | |
|---|---|---|---|---|
| edat | 1.093 | 1.049 | 1.139 | <0.001 |
| dbp | 1.062 | 1.012 | 1.114 | 0.014 |
| ecg: Anormal | 17.22 | 4.818 | 61.52 | <0.001 |
| ecg: Frontera | 3.466 | 1.021 | 11.76 | 0.046 |
> survminer::ggforest(mod) 
El millor model és el que conté les variables edat, dbp i ecg, ja que com hem vist en l’anàlisi univariant a edats i dbp elevades, i ecg anormal hi ha més probabilitat de morir abans.
3.2.2 Criteri clínic
El següent model s’ha ajustat per criteri purament clínic, és a dir, ajustem el model amb les variables desitjades. Com per exemple, volem saber si la supervivència entre els grups de ecg canvia si afegim al model la variable edatdiag.
> sel_var <- c("ecg", "edatdiag")
> frml <- formula(paste0("my.surv ~ ", paste0(sel_var, collapse = " + ") ))
>
> #N
> mod <- coxph(frml, data = dat)
> sum.surv <- summary(mod)
> c_index <- sum.surv$concordance
>
>
> print(kable(prettify(summary(mod), digits = 4)[,c(" ", "Hazard Ratio", "CI (lower)", "CI (upper)","Pr(>|z|)" )], booktabs = T,
+ caption = paste0( "Hazard Ratio Cox Multivariate" ,". Surv. C Index:",round(c_index["C"],2)), longtable = TRUE, escape = F, col.names = c(" ", "Hazard Ratio", "CI (lower)", "CI (upper)","P-value" )) %>%
+ kable_styling(latex_options = c("striped","hold_position", "repeat_header"), font_size = 14) %>%
+ row_spec(0,background = "#993489", color = "white") )| Hazard Ratio | CI (lower) | CI (upper) | P-value | |
|---|---|---|---|---|
| ecg: Anormal | 12.98 | 3.999 | 42.10 | <0.001 |
| ecg: Frontera | 2.308 | 0.7304 | 7.295 | 0.154 |
| edatdiag | 1.063 | 1.023 | 1.105 | 0.002 |
> ggforest(mod)